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ABSTRACT 

Aims. The simulation of three-wave interaction based plasma emission, thought to be the underlying mechanism for Type 111 solar 
radio bursts, is a challenging task requiring fully-kinetic, multi-dimensional models. This paper aims to resolve a contradiction in past 
attempts, whereby some studies indicate that no such processes occur. 

Methods. We self-consistently simulate three-waved based plasma emission through all stages hy using 2D, fully kinetic, electromag¬ 
netic particle-in-cell simulations of relaxing electron beams using the EPOCH2D code. 

Results. Here we present the results of two simulations; Run 1 (ni,lno = 0.0057, VhlAvi, = Vt/Ve = 16) and Run 2 {ni,/no = 0.05, 
Vfj/Avi, = Vij/Vi; = 8), which we find to permit and prohibit plasma emission respectively. We show that the possibility of plasma 
emission is contingent upon the frequency of the initial electrostatic waves generated hy the bump-in-tail instability, and that these 
waves may be prohibited from participating in the necessary three-wave interactions due to frequency conservation requirements. In 
resolving this apparent contradiction through a comprehensive analysis, in this paper we present the first self-consistent demonstration 
of fundamental and harmonic plasma emission from a single-beam system via fully kinetic numerical simulation. We caution against 
simulating astrophysical radio bursts using unrealistically dense beams (a common approach which reduces run time), as the resulting 
non-Langmiur characteristics of the initial wave modes significantly suppresses emission. Comparison of our results also indicates 
that, contrary to the suggestions of previous authors, an alternative plasma emission mechanism based on two counter-propagating 
beams is unnecessary in an astrophysical context. Finally, we also consider the action of the Weibel instability which generates 
an electromagnetic beam mode. As this provides a stronger contribution to electromagnetic energy than the emission, we stress 
that evidence of plasma emission in simulations must disentangle the two contributions and not simply interpret changes in total 
electromagnetic energy as evidence of plasma emission. 
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1. Introduction 

The emission of electromagnetic radiation at the local electron 
plasma frequency and twice the plasma frequency, otherwise 
known fundamental and harmonic plasma emission, has been 
observed to be a common phenomenon in both astrophysical 
and laboratory plasmas. Astrophysical examples include Type II 
& III solar radio bursts (Lin et al. 1981, 1986; Dulk et al. 1984; 
Goldman 1983; Robinson et al. 1993a,b, 1994; Reid & Ratcliffe 
2014), outer heliospheric radio emission (Kurth et al. 1984), 
and in planetary electron foreshocks (Etcheto & Eaucheux 
1984; Moses et al. 1984; Euselier et al. 1985; Lacombe et al. 
1985). Related processes in laboratory plasmas have also 
been discussed by, e.g. Hutchinson et al. (1978); Benford et al. 
(1980); Whelan & Stenzel (1981); Intrator et al. (1984); 
Whelan & Stenzel (1985). Eurthermore, radio bursts of mil¬ 
lisecond duration were recently discovered in the 1 GHz band 
Loeb et al. (2014), with strong evidence that they come from 
~ 1 Gpc distances, implying extraordinarily high-brightness 
temperature. Lyubarsky (2014) proposed that these bursts could 
be attributed to synchrotron maser emission from relativistic, 
magnetized shocks. Since fast radio bursts are associated with 
flaring stars (Loeb et al. 2014) it may well be that the plasma 
emission caused by an electron beam (as opposed to shocks) 
may be a viable mechanism. 


A long-standing, multiple-stage model for the fundamental 
and harmonic emission based on subsequent nonlinear three- 
wave interactions, has been considered extensively by authors 
including, but not limited to, e.g., Ginzburg & Zhelezniakov 
(1958); Melrose (1980, 1987); Cairns (1987); Robinson et al. 
(1994). In the first stage of this model, electron beams which are 
excited in the low-corona are susceptible to the bump-in-tail in¬ 
stability as they propagate through the background plasma, and 
thus the beams can generate Langmuir waves (L). The Langmuir 
waves can then decay into ion-sound waves (5) and electromag¬ 
netic waves at the local plasma frequency {T\) in the process 
L —> S + Ti, which is thought to explain fundamental emis¬ 
sion. Harmonic emission occurs due to two subsequent three- 
wave interactions; firstly, the electrostatic decay of forward- 
propagating Langmuir waves into ion-sound and backward- 
propagating Langmuir waves (L') L —> L' -H 5 (or equivalently, 
backscattering off ion-sound waves L + S ^ L'), and secondly, 
the coalescence of counter-propagating Langmuir waves to pro¬ 
duce electromagnetic waves at twice the local plasma frequency 
(T 2 ), L + L' ^ T 2 . 

Whilst many numerical works have considered individual 
stages of these fundamental and harmonic emission mecha¬ 
nisms, only a few attempts have been made to self-consistently 
verify the processes using fully-kinetic numerical simulations. 
The primary reason for this is that the self-consistent simula- 
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tion of the emission process is necessarily a computationally de¬ 
manding task that requires a full electromagnetic treatment in 
two dimensions which must be sufficiently large to contain all 
participating wavelengths. Two dimensions are required because 
the emission formulas (three-wave interaction probabilities) for 
both L ^ S + Ti and L + L' —> T 2 are only non-zero when 
a vector-product exists between participating wave vectors, 
see e.g., Eqs. (26.24) and (26.25) from Melrose & McPhedran 
(2005). As the instability time scales with the inverse of the 
beam-to-background density ratio (e.g., the quasilinear relax¬ 
ation timescale is Tqi = (nolni,) (vi,/Avi)^ direct simulation 
of astrophysical parameter regimes where beams are typically 
diffuse (riblnQ ss 10“^ - 10“®) requires the simulation to run for 
tens-to-hundreds of thousands of electron plasma periods. Fur¬ 
thermore, this is compounded in that huge particle numbers are 
required to correctly resolve the expected quasilinear relaxation 
timescales with the particle-in-cell (PIC) method (Ratcliffe et al. 
2014; Lotov et al. 2015). Lotov et al. (2015) found empirically 
that for beam-plasma systems the number of pseudoparticles per 
cell required scales as the inverse of the fraction of (real) par¬ 
ticles which are in resonance with the beam. In particular, for 
diffuse astrophysical beams, this leads to a very high pseudo¬ 
particle requirement of tens-to-hundreds of thousands particles 
per cell. 

In the handful of papers that have considered the fully self- 
consistent problem with the PIC method, a comparative review 
suggests that the previous results are somewhat contradictory. 
In the most extreme, some authors interpret the results of their 
experiments as evidence supporting the three-wave based mech¬ 
anisms (Kasabaetal. 2001; Umeda 2010), whereas others re¬ 
port that the such processes do not proceed in their simulations, 
and rather must invoke different mechanisms to produce emis¬ 
sion; such as the requirement for a counter-propagating beam 
(Ganse et al. 2012b,a, 2014; Timofeev & Annenkov 2014), or 
linear mode conversion off density inhomogeneities as per 
Sakai et al. (2005). Although, note that in the latter work they 
did not present results for a homogeneous setup, so arguably the 
underlying process responsible for emission in their work is not 
conclusively demonstrated. 

More subtly, further ambiguity exists in the details of the 
simulations proporting to show evidence of plasma emission. 
Firstly, limited attention is paid towards distinguishing the peaks 
in spectral power at the expected frequencies, which are taken 
as evidence of emission occurring, from thermal and noise lev¬ 
els. This is of crucial importance as noise organising itself as 
electromagnetic waves (most obviously manifest as power en¬ 
hancements along dispersion curves in (a>, k)-diagrams) is a typ¬ 
ical and unavoidable feature of PIC simulations. Secondly, pseu¬ 
doparticle numbers used are typically low (e.g., Kasaba et al. 
(2001) uses 16 particles per cell for background partxicles and 
4 for the beam particles, and Umeda (2010) uses 256 parti¬ 
cles per cell per species), and few papers report convergence 
testing of results. Umeda (2010) reported that increasing pseu¬ 
doparticle numbers diminished the signal associated with har¬ 
monic emission. Under-resolved particle numbers could con¬ 
tribute to excessive noise on the EM dispersion curves (and be 
mis-interpreted as emission as per the first point) and also result 
in poor replication of physical time-scales as per Ratcliffe et al. 
(2014) and Lotov et al. (2015). Finally, all of the previous works 
consider systems with high beam-to-background density ratios 
(all take Wi/no > 1%) and are so in the strong beam regime. 
As astrophysical systems typically have diffuse electron beams, 
this is presumably motivated by the faster relaxation times for 
dense beams allowing for significantly reduced computational 


Table 1. Simulation Parameters 


Run 

nblno 

Vb/Avb - VblVe 

Notes 

1 

0.0057 

16 

Weak beam 

2 

0.05 

8 

Dense beam 

3 

0.05 

8 

Two counter- 
propagating beams 


requirements. However, substantial theoretical work has repeat¬ 
edly shown that beam-plasma systems are expected to posses in¬ 
credibly sensitive parameter spaces. In particular, dense beams 
are associated with strong modification of the dispersion rela¬ 
tionships (O’Neil & Malmberg 1968; Cairns 1989), and so we 
raise the question how valid are strong beam simulations of a 
process that is analytically (in quasilinear theory) prescribed for 
weak beams? 

The primary aim of this paper is to carry out detailed 2D PIC 
simulations of varying beam-plasma systems to resolve the ap¬ 
parent contradictions discussed above; in other words to clarify 
whether the three-wave interaction based, single-beam plasma 
emission can be self-consistently shown to proceed as expected 
for astrophysical parameter regimes, or whether the two-beam 
mechanism is required as suggested by Ganse et al. (2012b). In 
resolving this apparent contradiction through a comprehensive 
analysis, in this paper we present the first self-consistent demon¬ 
stration of fundamental and harmonic plasma emission from a 
single-beam system via fully kinetic numerical simulation. The 
paper is structured as follows; in Section 2 we detail the numer¬ 
ical methods used and the initial conditions/setup for the exper¬ 
iments, in Section 3 we present the results, in Section 4 we dis¬ 
cuss their implications and in Section 5 we draw concluding re¬ 
marks. 

2. Numerical Setup 

The simulations are carried out using EPOCH2D, a 2.5D Bird- 
sail and Langdon type PIC code with a core solver based on the 
PSC code of Rhul (Chapter 2 of Bonitz & Semkat (2006); see 
also Brady & Arber (2011)). In this paper we present the results 
of three numerical experiments. The first two simulations vary 
in beam-to-background number density ratio, beam temperature 
and speed as Wi/no = 0.0057, Vi/Av* = VblVg = 16 (Run 7) 
and nblno - 0.05, Vbl?svb - VblVe - 8 (Run 2). Thus, Run 1 
is weaker beam than considered in previous PIC simulations of 
plasma emission (but still orders of magnitude larger than typi¬ 
cal astrophysical plasmas) and Run 2 is more comparable to the 
past simulations. The initial conditions are chosen for two rea¬ 
sons; firstly, their ID analogues have recently been considered 
by Baumgartel (2014), who noted a difference in the nature of 
the electrostatic decay processes in the two regimes, and so we 
may explore the consequences for emission processes in our ex¬ 
tension to 2D. Secondly, each set-up gives comparable parame¬ 
ters 

P = (n,/no)‘^\vi/Av,) (1) 

of P = 2.85 and F = 2.94, which is a measure of the kinetic ver¬ 
sus fluid nature of the instability. As per the parameter study of 
Cairns (1989), these F values indicate a predominantly kinetic 
instability (NB: a key finding of that paper was that the kinetic 
criterion P < 1 sometimes reported by other authors is not accu¬ 
rate). These experiments are presented in order to illustrate the 
sensitivity to different beam parameters, and to investigate the 
consequences for efficiency of plasma emission. Furthermore, a 
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Fig. 1. Electron velocity distribution functions for Run 1 (a-c) and Run 2 (d-f) at times = 0, 300, 600 respectively. The solid line shows the 
velocity distribution aligned with the beam {vx) and the dashed line shows the transverse velocity (v^,). Panels (g) and (h) are representative of the 
behaviour in 2D phase-space for Run 1 and 2 respectively at = 300. Note the plateau formation in with the saturation of the bump-in-tail 
instability, and the transverse heating of the beam electrons due to the Weibel instability. 


third experiment {Run 3) extends Run 2 by considering the setup 
in the presence of an additional counter-streaming (but otherwise 
identical) electron beam in order to investigate the differences 
between single-beam and two-beam plasma emission, and the 
circumstances in which a two-beam emission mechanism may 
be necessary. The difference in parameters between the three 
runs is summarised in Table 1 

All three experiments share common background plasma pa¬ 
rameters, which are chosen as follows; the background number 
density no = «« = «!' = 4 x 10® m“^, background temperature 
Te - 200 eV, Ti = 0J3Te and the mass ratio niilnie - 1836. The 
background magnetic field is set to zero (i.e., we do not consider 
a weak guide field aligned with the beam). The simulations are 
ran for lOOOw^] in a domain of size - Ly - 600Tz), where 
the cell width is equal to the Debye length - t^y - Ao, and 
periodic boundary conditions are applied. For each simulation, 
convergence testing was carried out for pseudoparticle counts of 
125, 250, 500 and 1000 particles per cell per species (PPCPS). 


The convergence tests consisted of (1) ensuring convergence of 
relaxation time of the bump-in-tail instability (i.e., ensuring con¬ 
vergence for figures similar to Fig. 1) and (2) ensuring conver¬ 
gence of spectral energy density levels associated with different 
wave modes to within an order of magnitude (NB: convergence 
to specific numbers is impossible due to the random nature of the 
PIC code; increasing the particle count improves statistical rep¬ 
resentation). The convergence tests found tolerable convergence 
of the results for 500 and 1000 PPCPS. The results presented 
here are for 1000 PPCPS, which required ~ 12 hours of comput¬ 
ing time each on 480 2.50GHz Intel Xeon cores. 


3. Results 

First we consider the evolution of the electron population in ve¬ 
locity phase space. Figure 1 shows electron velocity distribu¬ 
tion functions for Run 1 (a-c) and Run 2 (d-f) at times ttOpe - 
0, 300, 600 respectively, where the solid line shows the veloc- 
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Fig. 2. Fast fourier transforms of for Run 1 (left) and Run 2 (right) over times 50 < < 150, 450 < < 550, 850 < tai^j, 

top to bottom respectively. 
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Fig. 3. Spectral energy density (J m along the Langmuir wave disper¬ 
sion curve for Run 1 (top) and Run 2 (bottom) for subsequent 100 
time windows (centred ±50(i>^l about the indicated time). Note the de¬ 
velopment of a broad spectrum in negative in Run 1 with energy 
densities several orders of magntude larger than the negative spec¬ 
trum of Run 2. The dotted lines show the expected initial value of the 
initially growing electrostatic modes. 


ity distribution aligned with the beam {v^) and the dashed line 
shows the velocity distribution across the beam (v,.). In f{v^) 
phase space (solid lines), we see the saturation of the initial 
bump-in-tail instability and the merging of the beam electrons 
to the bulk distribution, giving the characteristic plateau forma¬ 
tion. This is observed to occur by around f =» 60 (Run 1) and 
t X 600 ojpl (Run 2), which are in excess of time-scales predicted 
by the quasi-linear formula. We however, confirmed that this is 
a realistic relaxation time in our convergence testing, rather than 
being due to poor particle counts as reported by Ratcliffe et al. 
(2014). The discrepancy is due to the use of the relatively strong 
beams (in terms of kinetic energy); the applicability of the weak 
turbulence/quasi-linear results require that the ratio of beam to 
background kinetic energies be much less than unity. In f{vy) 
phase space (dashed lines) we see the heating of the electrons in 
both runs, due to the action of the Weibel instability. The Weibel 
instability is expected to cause perpendicular heating in beam- 
plasma systems except in the case of strong magnetic helds or in 
the reduction to a ID system, as discussed in detail by Karlicky 
(2009). Again, this occurs on a faster timescale and saturates 
faster in Run 2 compared to Run 1. Panels (g) and (h) of Fig¬ 
ure 1 show visualisations of the above behaviour in 2D phase- 
space at a representative time(to^j = 300) for Run 1 and Run 2 
respectively. These panels highlight that it is the beam popula- 



Fig. 4. kr Spectra of ion density fluctuations (integrated in y and t) for 
Run 1 (solid line) and Run 2 (dashed line). Note the presence of peaks 
at |k| Ikt, where kt is the predominant wavenumber associated with 
the growing Langmuir waves found in in Run 1, suggestive of the action 
of Langmuir decay / scattering processes involving ion-acoustic waves. 
No analogous peak for Run 2 is found. 

tion, rather than the background electrons, which are susceptible 
to heating in the y-direction, which has consequences that we 
will explore later. Overall, for both runs, we see the relaxation 
of electron phase space to asymptotic states during the course of 
the experiment. 

We now turn our attention to electrostatic waves generated 
by the beam relaxation. Figure 2 shows the evolution of the par¬ 
allel electric held energy in (kx,a>) space by considering 2D, 
windowed Fourier transforms in [x, t) space over subsequent 
100 Wpp periods, which are then integrated over the y direc¬ 
tion. The time window size was determined by experimenta¬ 
tion and provides the best balance between frequency resolu¬ 
tion (Aw) and time cadence, allowing us to track the evolu¬ 
tion in both Fourier space and time, whilst preserving suffi¬ 
cient frequency resolution to reasonably compare spectra to ex¬ 
pected dispersion curves. The left-hand column shows the spec¬ 
tral electrical energy density for Run 1, and the right for Run 
2, over time periods 50 < tojpe < 150, 450 < tcjpe < 550 and 
850 < tcjpe < 950 (top to bottom). The dispersion curves for the 
expected beam modes (cj/k = Vh) and (unmodified) Langmuir 
modes (o/ = + 3k^Vg) are overlaid. 

For Run 1, we And that the majority of initial growth oc¬ 
curs on the intersection of the beam-mode and Langmuir mode. 
Over time, the power is efficiently transferred to negative k^ via 
decay/backscattering processes of the type L ^ L' + S , thus 
creating a seed population of forward and backward propagat¬ 
ing waves which are possible candidates for participation in the 
coalescence L + L' —> T 2 required for harmonic emission. The 
backwards (negative k^) portion of the Langmuir wave disper¬ 
sion curve appears to be slightly shifted towards higher oj at 
larger k^ than expected in a quiescent plasma, which is consistent 
with the solutions of Cairns (1989), who numerically determined 
the changes to plasma dispersion relation due to the presence of 
a beam. We can more closely inspect and compare the energy 
evolution associated with electrostatic waves by plotting the evo¬ 
lution of energy density along the curve u/ = + 3k^V^, as 
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shown in Fig. 3, which illustrates the previous discussion. Note 
that, in Fig. 3, we have sampled the frequency range +Aw above 
and below the curve to account for observed deviations in co 
from the unmodified Langmuir wave. The role of three-wave 
backscatter/decay progress of the type L —> L' -H 5 in the cre¬ 
ation of the counter-propagating Langmuir wave propagation is 
confirmed by the presence of enhanced ion-density fluctuations 
at the expected ks (ks ~ 2kL, as shown in Figure 4). Further, 
this is consistent with the PIC simulations of Baumgartel (2014) 
who considered the ID analogue of this system (i.e., the same 
beam-plasma setup as Run 1 in a ID PIC code) which afforded 
superior kx resolution by using box sizes which are inaccessible 
to our 2D simulations, unambiguously confirming the action of 
decay/backscattering processes in the the system’s ID counter¬ 
part. 

Thus, for Run 1, we find qualitative evolution of the power 
spectra that is consistent with the first two stages of the single¬ 
beam, three-wave based harmonic emission mechanism, namely 
(1) the efficient coupling of the beam to electrostatic waves, and 
subsequently (2) the generation of an eligible seed population of 
forward and backward propagating electrostatic waves via non¬ 
linear three-wave interactions. Quantitatively, the seed popula¬ 
tion grows to be at least 2-3 orders of magnitude in excess of 
background levels by tojpe = 900. Additional features can also 
be identified, most notably the so-called high-harmonic electro¬ 
static nonlinear plasma waves (Rhee et al. 2009; Yi et al. 2007; 
Yoon et al. 2003). As pointed out by Yoon et al. (2003) these 
modes are nonlinear eignemodes that exist by virtue of the non¬ 
linear plasma response to the presence of the beam, and do not 
exist as ‘natural’ modes in quiescent plasmas. We also note the 
presence of a spectral peak at kx - 0. This corresponds to a 
beam-aligned, standing mode of the electric field (Ex) oscillating 
at the local plasma frequency, which is present from the simula¬ 
tion initialisation (thus, before the instability onset). It is caused 
by the non-zero initial current imposed by the beam at f = 0, and 
has been discussed in detail by Baumgartel (2013). Whilst it is 
possible to remove this mode by introducing compensating drift 
velocity to the background electrons, we tolerate its presence as 
it is unclear whether such a compensation is physically appro¬ 
priate (in particular, it may influence the correct return-current 
processes). Regardless, we have found that the amplitude asso¬ 
ciated with this mode, in all runs was much less than that of the 
Langmuir waves generated after the instability onset and as such 
it does not go on to effect the dynamics of the Langmuir waves 
which participate in the plasma emission mechanism. The influ¬ 
ence and relative amplitude of the beam-aligned mode is most 
easily determined from time-distance diagrams (not shown here, 
for brevity), however its minor influence on the energy budget 
can be seen Figure 8, which is discussed in more detail later in 
this section. 

However, in Run 2 we do not find analogous behaviour. The 
beam mode ‘decouples’ from the expected Langmuir / elec¬ 
trostatic mode and the main concentration of power associated 
with the beam mode drifts towards higher kx as time evolves 
(Fig. 2). Furthermore, it is primarily concentrated below the 
frequency expected of forward-propagating electrostatic waves, 
and so there exists the possibility of three-wave interactions be¬ 
coming strongly limited by frequency conservation requirements 
as shown by Cairns (1989); this is discussed further in Section 4. 
Growth on both the forward and backward portions of the Lang¬ 
muir curve itself (Fig. 3) is predominantly at small \kx\. At all 
times, the backward (negative kx) spectrum is at least an order 
of magnitude weaker than possible counterparts in Run 1 (but 
generally ranges from being 1-3 orders weaker), supporting 


the conjecture that the decay process has been inhibited in the 
case of Run 2. Additionally, an absence of significant ion-density 
power enhancement in Run 4 is consistent with an absence of 
decay/scattering processes in excess of noise levels, contrary to 
Run 1 (cf. Fig 4). 

We now turn our attention to the consequences of the no¬ 
tably different nature of electrostatic wave growth and evolution 
in Runs 1 and 2 for plasma emission by considering equivalent 
figures for the transverse magnetic field B^. In Figure 5 we see 
the evolution of spectral energy density in in Run 1 (left) 
and Run 2 (right) over the same timesteps considered in Fig¬ 
ure 2. For both runs we see the growth of power along the elec¬ 
tromagnetic dispersion curves (overplotted), some of which is 
associated with the PIC noise discussed in Section 1. When in¬ 
vestigating plasma emission processes with the PIC method, it 
is crucial that any signals of emission at the expected wavenum¬ 
bers and frequencies is clearly distinguished from such noise. 
For the results presented here, we estimate the noise threshold 
from the average power on the curves away from the frequencies 
that are expected to be enhanced by plasma emission mecha¬ 
nisms (a> - cOpe, and its harmonics), and find it to be of the order 

-140dB = 10“^^^ Jm“^. Hence forth, we only consider signal 

on electromagnetic dispersion curves that is distinguished above 
this threshold as an indicator of plasma emission occurring. 

With regards to indicators of fundamental emission in Fig¬ 
ure 5, in both experiments we see an enhancement of power at 
to - tOpe for small kx- We can take a closer look at the enhance¬ 
ment in Figure 6, which considers the distribution of spectral en¬ 
ergy density over kx?Ato - cOpe (similar to Figure 3). At small kx 
for Run 1, we observe the growth of a peak in excess of the noise 
levels, which reaches ~ 10 '^Jm“^, two orders of magnitude 
above the threshold, in the time window 900 < ttOpe < 1000. 
The peak’s apex is slightly shifted towards negative kx, located 
at kxAo - -0.01 . Note that, the expected wavenumber from the 
fundamental emission process L —> 5 -H Ti for Run 1 ’s beam 
parameters is kriAo ~ 0.002, below our sampling rate of A^. 
As such, confirming whether the signal is predominantly con¬ 
centrated at the expected wavenumber is not possible for these 
simulations and would require larger spatial domains. For Run 
2, we see the development of a weaker, narrower signal peaked 
at kx = 0 which is at most an order of magnitude above the 
noise threshold, i.e., its peak spectral energy density is one or¬ 
der of magnitude below the equivalent signal in Run 1. Thus, 
for both runs we see a signal that is consistent with fundamen¬ 
tal emission, but it is comparatively weaker for Run 2. Addi¬ 
tionally, for both runs, in Figures 5-7 we see features which 
are not associated with electromagnetic radiation (i.e, the peak 
at larger positive kx in figure 6). These are driven by the per¬ 
pendicular heating of the beam electrons (Weibel instability) as 
seen in Figure 1 and are manifest its various harmonics (cf. Fig¬ 
ures 2 and 5). These electromagnetic beam modes are found in 
the Ey spectrum (which is not shown here, but is qualitatively 
as per Fig. 5). On close inspection, these modes identified can 
be seen in other authors works (eg., Plate 2 of Kasaba et al. 
(2001)) although they recieved no attention at the time, pre¬ 
sumably due to negligence in considering perpendicular phase 
space, and so the realisation of the role of the Weibel insta¬ 
bility, which was only first demonstrated by Karlicky (2009). 
The tenancy for the Weibel instability to drive non-maxwellian, 
electromagnetic eignemodes at to - kvb was discovered experi¬ 
mentally by Urrutia & Stenzel (1984) and explored theoretically 
by Goldman & Newman (1987). To our knowledge, this is the 
first time the the electromagnetic beam mode has been discussed 
in the context of plasma emission. Relatively large amounts of 
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power are concentrated in these modes (more than that emitted 
by the plasma emission) and so we caution that total power con¬ 
tained in the electromagnetic field is not exclusively associated 
with the plasma radio emission and must not be interpreted as 
such. This will be important in our later consideration of the 
overall energy budget of the systems. 

In Figure 5, as time advances we see the development of a 
clear enhancement near u) - 2(Ope on the electromagnetic disper¬ 
sion curve for Run 1, but equivalent signal for Run 2 can only be 
poorly distinguished against the EM noise levels. Again, this is 
examined more closely in Figure 7 by considering the spectra at 
fixed CO - 2(0 pe. For Run 1, we find four clearly distinguished 
spectral peaks close to the expected \kx\ values for harmonic 
emission generated by the coalescence L + L' —> T 2 which grow 
to be in excess of two orders of magnitude above noise levels 
(ss 10“^^ Jm“^) by the time period 900 < t(Ope < 1000. For Run 
2, however, we only see the slow development of small peaks 
that are comparable with the noise, with a maximum spectral 
energy density of ss 10“^"^Jm“^. Thus, signals consistent with 
harmonic plasma emission that are clearly distinguished from 
background levels are detectable in only the case of Run 1. Note 
that the reason that we see four peaks in at w = 2cOpe, rather 
than two, is simply an artefact of projecting obliquely propaga¬ 
tion emission into kx space. If we instead consider B^(kx, ky) at 
fixed CO = 2cOpe, we find power enhancements suggesting emis¬ 
sion at angles (relative to beam direction) of 22°, 68°, 112°, and 
158° which is in agreement with that reported by Umeda (2010). 

In Figure 8 we consider the energy evolution of the sys¬ 
tem, normalised to the initial beam kinetic energy. For the 
energy associated with the parallel electrostatic component 
(f O.SsqE^ dV), for both runs, we initially see the ‘beam-aligned 
mode’ (kx - 0 mode). After the instability onset, the energy as¬ 
sociated with Ex grows to around 10% of the beam kinetic en¬ 
ergy in both cases. Thus, the driven electrostatic energy rapidly 
overpowers the kx -0 mode, as discussed earlier. This is accom¬ 
panied by an energy increase associated with Ey (J O.SsoEy dV), 
which is predominantly by the perpendicular acceleration of the 
beam particles by the Weibel instability, but will also contain 
the electric energy associated with any plasma emission. The 
energy associated with the corresponding magnetic field com¬ 
ponent, (J 0.5/to'Bj dV), similarly grows after the instabil¬ 
ity onset, then decays, and subsequently begins to grow again at 
later times in both runs. Here we stress that the initial growth 
is associated with the direct generation of electromagnetic beam 
modes by the Weibel instability. The later growth is associated 
with the comparatively weak fundamental and harmonic emis¬ 
sion processes identified in earlier paragraphs. Here, we stress 
caution against interpreting any increase in total electromag¬ 
netic energy as evidence of emission as has been the case in 
some previous studies. For example, from Figure 8, one may 
be tempted to conclude a similar efficiency and timescale of 
presumed ‘emission processes’ for both runs, although as we 
have seen in our earlier analysis that emission, particularly har¬ 
monic, is much weaker or non-existent in Run 2. It is crucial to 
appreciate that the electromagnetic beam mode is the predom¬ 
inant contributor to total electromagnetic energy, and so emis¬ 
sion energy and efficiency is better identified by the spectral 
methods we have employed which are able to disentangle the 
contributions of different modes. We also note that energy in 
Bx, By and does not grow beyond noise levels and as such 
contain only a negligible portion of the systems energy. From 
this we may infer that emission is linearly polarised, as it is 
only manifest in and not By (equivalently, manifest in Ey and 


not Fj.). This is to be expected, due to the absence of an ini¬ 
tial magnetic field. Plasma emission from Type III solar radio 
bursts has be observed to be typically weakly circularly polarised 
(McFean 1971; Suzuki & Sheridan 1977; Dulk & Suzuki 1980; 
Suzuki & Dulk 1985), with harmonic emission having a lower 
degree of polarisation than fundamental. Finally, we note that 
when a total system energy is calculated we find that it is well- 
conserved, to an accuracy of 0.03% during the simulation life¬ 
time. 

4. Discussion 

We have found the following key features in the comparison of 
the wave dynamics and consequences for plasma emission in two 
different (single) beam-plasma systems; 

1. The growth of the beam mode on the Fangmuir wave disper¬ 
sion curve in Run 1, and apparent ‘decoupling’ of the beam 
and Fangmuir modes in Run 2, whereby the presence of the 
denser beam has significantly modified the nature of the for¬ 
ward propagating electrostatic wave modes. 

2. The apparent difference in susceptibility of the two systems 
to processes of the type L ^ L' + S, with consequences for 
the production of a seed population of counter-propagating 
Fangmuir/Electrostatic waves. 

3. We observe fundamental emission distinguishable above the 
noise threshold in both cases. At its peak, the spectral energy 
density associated with the emission is two-orders of mag¬ 
nitude stronger than the noise threshold in Run 1, but only 
one-order of magnitude greater in Run 2. Thus, the funda¬ 
mental emission is weaker in the strong beam case of Run 
2 . 

4. We observe that harmonic emission is clearly distinguishable 
above the noise level for Run 1 (two orders of magnitude at 
peaks), but only weak signals comparable to background lev¬ 
els are detected in Run 2. Thus, harmonic emission is much 
less effective in the case of Run 2, possibly to the point of 
being prohibited altogether. 

That statement (4) is the case, namely that harmonic emis¬ 
sion is only observed at twice the background level in Run 2, 
and that it is in its most generous interpenetration weak in Run 
2, is a direct result of (2). Run 1 has been demonstrated to posses 
a larger seed population of forward and backward propagating 
electrostatic waves for participation in the processes of the type 
L + L' —> 7’2. Thus, as is it may be expected, for harmonic emis¬ 
sion the efficiency of electrostatic decay processes L L' + S 
is of primary importance. The fundamental difference in the ef¬ 
ficiency of the transfer of energy from the beam mode to a pop¬ 
ulation of forward and backward propagating Fangmuir waves 
can be explained in terms of (1), i.e. the nature of the gener¬ 
ated forward propagating electrostatic waves in the two parame¬ 
ter regimes. Crucially, in Run 2 the forward propagating electro¬ 
static mode associated with the beam is primarily concentrated 
below the plasma frequency (cf. Fig. 2). Reading from Figure 2, 
the packet is concentrated at approximately 0.9 cOpe- To partici¬ 
pate in the three wave interactions L L' ±S (decay/scattering) 
we require 0.9 ^ cou IcOpe ± cos jcope, where cos is a frequency 
available to ion-sound waves and cou is a frequency available to 
the backwards branch of the Fangmuir / electrostatic waves. For 
simplicity we can simply note the requirement coulcope > 1, 
and so we require the presence of ion-sound waves with fre¬ 
quencies cosjcope >0.1, which is beyond the cut-off frequency 
(co pi!COpe ~ 0.02). Thus, for Run 2, such interactions are strongly 
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Fig. 5. Fast fourier transforms of for Run 1 (left) and Run 2 (right) over times 50 < < 150, 450 < < 550, 850 < < 950, from 

top to bottom respectively. 


Article number, page 8 of 12 















































Thurgood & Tsiklauri: PIC Simulations of Plasma Emission 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

kAo 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

kAo 


Fig. 6. Spectral energy density in B, (Jm“^) at the fundamental fre¬ 
quency (jj = LL)pg for Run 1 (top) and Run 2 (bottom) for subsequent 
time windows. 

suppressed as most of the beam mode occupies prohibited fre¬ 
quencies. On the other hand, for Run 1 the majority of the for¬ 
ward propagating waves are concentrated on the Langmuir dis¬ 
persion curve and so they can interact with the low-frequency 
ion-sound waves, resulting in the development of the broad spec¬ 
trum of counter-propagating waves as per (2), which go on to co¬ 
alesce into harmonic emission as per (4). Our demonstration of 
this decay inhibition due to breaking of frequency matching re¬ 
quirements (wave beat conditions) is the first verification of the 
arguments made from a theoretical perspective by Cairns (1989). 

The observed reduction in fundamental emission efficiency 
as per (3) can be explained similarly. Fundamental emission pro¬ 
cesses based on three-wave interactions of the type L± S Ti 
demand frequency be conserved such that oji ± a>s = cot, ■ As 
the beam mode lies signihcantly below the a>pe, electromagnetic 
waves must be above it, the frequency conservation requirements 
require the low-frequency wave to be in excess of the cut-off, and 
thus they cannot be met. 

In summary, we have presented two experiments for appar¬ 
ently comparable beams (at least in terms of the reactive/kinetic 
P parameter, and a difference in nt/no of only one order of mag¬ 
nitude), whereby only in the case of Run 1 do we see significant 
and unambiguous fundamental and harmonic plasma emission 
far in excess of background levels. We propose that the demon¬ 
strated sensitivity, particularly to beam density, is the underly¬ 
ing explanation as to why Ganse et al. (2012b) have reported be¬ 
ing able to find no evidence consistent with plasma emission in 
PIC simulations. Furthermore, we caution that as previous stud¬ 



kAo 
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Fig. 7. Spectral energy density in (J m“^) at the harmonic frequency 
10 = IcOpi; for Run 1 (top) and Run 2 (bottom) for subsequent time 
windows. 

ies such as Kasaba et al. (2001); Umeda (2010) have interpreted 
their results as signs of emission in similar parameter regimes 
to Run 2 but not reported on noise thresholds, it is unclear if 
such systems actually generate efficient plasma emission that is 
compatible with observed radio emission. Indeed, typical setups 
(with dense beams) are closer to those used in Run 2 (i.e., the 
case where we find limited evidence for fundamental emission, 
and no evidence for harmonic emission.) 

Ganse et al. (2012b), who could not demonstrate three-wave 
based emission from a single beam, explored different mech¬ 
anisms to produce (harmonic) plasma emission. One notable 
mechanism, based on the action of two counter-propagating 
beams, was discussed in a series of papers (Ganse et al. 2012b,a, 
2014). To illustrate some key differences in emission arising 
from the (single-beam) fundamental and harmonic plasma emis¬ 
sion mechanisms and the two-beam emission mechanism, we re¬ 
ran Run 2 and included an additional, but otherwise identical 
counter-propagating electron beam (Run 3). 

From early times we see strong signals in that are consis¬ 
tent with harmonic emission (Fig. 9). This is explained in that 
the beams directly drive electrostatic modes that match the beat 
conditions for harmonic emission (Fig. 10). As such, the gener¬ 
ation of harmonic emission is guaranteed if the two beams have 
identical properties (as this determines their a> and |k|). It is un¬ 
clear how robust this mechanism is in the case of non-identical 
beams of different nt/no, Vi/Av/,, Vbjvp, and would make for in¬ 
teresting future work. We stress that this is not the ‘classical’ har¬ 
monic emission which relies on three subsequent stages (growth 
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Fig. 8. Energies normalised to initial beam kinetic energy for Run 1 
(top) and Run 2 (bottom). The black lines show the kinetic energy as¬ 
sociated with the beam and background electrons (in the x-direction), 
the red lines the electric field energy (from x, y, and z components when 
reading from top to bottom) and blue the magnetic field energy com¬ 
ponent at ~ 10“"*, and 6^ and By components at noise level ~ 10“’). 


of forward propagating electrostatic waves, backscatter and de¬ 
cay to produce a counter propagating population, and their coa¬ 
lescence), and is rather directly powered by the artificial initial 
conditions. 

We also note that, compared to Run 2, we find little evidence 
for enhanced fundamental emission owing to the presence of the 
additional beam. This is consistent with the argument presented 
as to why L —» S + Ti does not proceed efficiently Run 2; we 
observe in Fig 10 that the two beam modes are generated at the 
similar a> and k as in Run 2 and so we do not expect a three wave 
interaction which is prohibited in Run 2 to be permitted due to 
the additional beam. Finally, we comment on the low-frequency 
enhancements in the electrostatic power spectrum. These are not 
ion-sound waves as their frequency is too high, being above the 
cut-off frequency. They are rather a direct, nonlinear plasma re¬ 
sponse to the ponderomotive force of the two crossing beam 
modes. Apparent harmonics of these ‘daughter waves’ (which 
are actually driven by the corresponding harmonics of the beam 
mode), can also be observed in Fig 10. 


5. Conclusion 

We have presented two numerical experiments for different 
beam-plasma systems which demonstrate a remarkable sensi¬ 
tivity in terms of resulting wave dynamics, with drastic conse¬ 
quences for plasma radio emission. 



In Run 1 (a more tenuous and fast beam of n/j/no = 0.0057, 
Vfe/Avfc = VbjVe - 16) , we find evidence consistent with all 
stages of the three-wave based fundamental and harmonic emis¬ 
sion mechanisms, including the beam-mode to Langmuir mode 
coupling, the growth of a population of counter-propagating 
Langmuir/electrostatic waves via backscattering and decay pro¬ 
cesses, the action of fundamental emission, and the coales¬ 
cence of the counter-propagating population to produce har¬ 
monic emission. For the first time using a fully kinetic and elec¬ 
tromagnetic PIC simulation, we can confirm the role of all such 
stages whilst taking care to distinguish signals above the inher¬ 
ent noise levels associated with particle methods, and perform¬ 
ing appropriate convergence testing. Thus, Run 1 is arguably the 
first unambiguous confirmation of the three-wave based emis- 
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Fig. 10. Fourier transforms of for Run 3 at early (top, 50 < < 

150) and late (bottom, 850 < < 950) times. 
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sion processes resulting from a single electron beam in the liter¬ 
ature using the fully kinetic PIC approach. 

For Run 2 (a slower, denser beam of nt,lno - 0.05, VblAvb = 
Vb/Ve - 8), we demonstrate the sensitivity of the parameter 
space by considering a more dense beam with rii/no, a similar 
density to past works, and found that the processes are signif¬ 
icantly suppressed due to the resulting non-Langmuir charac¬ 
teristics of the beam mode. Whilst a full parameter study may 
prove useful, but is beyond available computational resources, 
we hope that Run 2 demonstrates plainly that caution must be 
applied when attempting to simulate astrophysical beam-plasma 
systems using unrealistically dense beams. Whilst a larger den¬ 
sity ratio reduces relaxation times (i.e. computer time), the re¬ 


sults are unlikely to be physically representative of the intended 
system due to the sensitivity to beam parameters. 

We also make the first connection between the action of 
the Weibel instability and the generation of an electromagnetic 
beam mode in the context of plasma emission. As this provides 
a stronger contribution to electromagnetic energy than the emis¬ 
sion, we stress that evidence of emission in simulations must dis¬ 
entangle the two contributions (such as by our spectral approach) 
and not simply interpret changes in total electromagnetic energy 
as emission. Following Karlicky (2009), who found that only 
very strong fields (cjpelojce ~ 1) could inhibit this transverse be¬ 
haviour, we expect that effect this is of importance to application 
to solar radio bursts (where the field is relatively weak). 

Finally, comparison of our results indicate that, contrary 
to the suggestions of authors including Ganse et al. (2012b,a, 
2014), the two beam, or counter-propagating beam, mechanism 
is not necessary to produce harmonic emission, and that in cer¬ 
tain parameter spaces (such as Run 1), the single-beam emission 
mechanisms can proceed. However, that is not to say that the 
counter-propagating beam mechanism does not work. Rather, 
the suitability of either process depends on the physical situation 
and crucially, whether two beams are expected. In cases where 
we expect two counter-propagating beams to exist which also 
happen to be suitably connected to the Langmuir mode (primar¬ 
ily, this implies lower density ratios than those considered by 
(Ganse et al. 2012b,a, 2014), we anticipate competition between 
the two different mechanisms. 
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